1 データの用意

今回は前々回と同じデータを生成して分析していく。

個人の所得\(Y\)と学歴\(X\)・能力\(A\)との関係についての架空データを次のように生成する（ただし\(\varepsilon\)は測定誤差などを反映した撹乱項）

\[ Y = 200 + 10A + 500X+ \varepsilon \]

  • サンプルは10万人
  • 切片は200万円
  • 能力が1上がると所得は10上昇する
  • 能力は0から100まで均等に分布する
  • 大卒だと所得が500万円上昇する
  • 大卒になる条件は2つあり，どちらかの条件を満たせば大卒になるとする
    • 条件1：能力が 80 以上の約 2 万人の中から約 1 万人がランダムに選ばれて大卒となる。
    • 条件2：能力を部分的に反映した学力テストの点数が 180 点以上であれば大卒となる（約1万人が合格する）
# 事前準備 --------------------
# パッケージの読み込み
library(tidyverse)
# 乱数の種を固定
set.seed(0)

# データの生成 ----------------
n = 100000
# 能力は0から100まで均等に分布
ability = runif(n, min = 0, max = 100)

# IDとabilityをデータフレームに格納する
df = data_frame(ID = 1:n, ability)

# 大卒フラグ
## 条件１：能力が 80 以上の約 2 万人の中から約 1 万人がランダムに選ばれて大卒となる。
college = df %>% filter(ability >= 80) %>% sample_frac(0.5) # 大卒の人
college["is_college1"] = 1
no_college = anti_join(df, college, by = c("ID","ability")) # 大卒じゃない人
no_college["is_college1"] = 0
df = bind_rows(college, no_college) %>% arrange(ID) # 両者を結合

## 条件２：能力を部分的に反映した学力テストの点数が 180 点以上であれば大卒となる
### １万人くらいが該当するようにする
df["score"] = 30 * log10(ability) + rnorm(n, mean = 115, sd = 10)
df["is_college2"] = 1*(df["score"] >= 180)

## ２つの条件のいずれかに該当していれば大卒とする
df = df %>% mutate(is_college = case_when(is_college1 == 1 | is_college2 == 1 ~ 1, # "|"はorの記号
                                          TRUE ~ 0)) # それ以外は0
# 所得
df["income"] = 200 + 10*df["ability"] + 500*df["is_college"] + rnorm(n, sd = 50)

# 最初の6行
head(df)
## # A tibble: 6 x 7
##      ID ability is_college1 score is_college2 is_college income
##   <int>   <dbl>       <dbl> <dbl>       <dbl>      <dbl>  <dbl>
## 1     1    89.7           1  172.           0          1  1572.
## 2     2    26.6           0  143.           0          0   479.
## 3     3    37.2           0  158.           0          0   492.
## 4     4    57.3           0  172.           0          0   720.
## 5     5    90.8           1  184.           1          1  1738.
## 6     6    20.2           0  170.           0          0   328.

こうして生成したデータの所得と能力・学歴・テストの点数の関係をプロットすると次のようになる。

以下では「学歴（大卒になること）が所得をどれだけ上昇させるのか」を知ることを目的として分析していくことにする。

2 回帰不連続デザイン

2.1 概要

temp <- df %>% 
         mutate(education = case_when(is_college2 == 1 ~ "大卒(条件2)",
                                      is_college == 0 ~ "非大卒")) %>% 
         filter(!is.na(education))

ggplot(temp,
       aes(x = score, y = income, color = education, fill = education))+
  geom_point(alpha = 0.5)+
  geom_smooth(method="lm",se=F, color="black", alpha = 0.5)+
  labs(title = "大卒（条件2）と非大卒のスコアと所得")

条件2で大卒になった人の場合を見ていくと，テストの点数がある閾値（180点）を超える点で所得にジャンプがみられる。 回帰不連続デザイン(regression discontinuity design: RDD)はこうした閾値を伴うDGPを利用する調査設計（デザイン）で，以下のようなロジックをもとに分析を行う。

  • スコアが180点になる点で大卒になれるかどうかを決めるというのは人為的な制度であり，その背景で個々人の能力は連続的に分布している（179点で不合格だった人と180点で合格した人とで能力に大差はない）と考えられる
  • 能力\(A\)が所得\(Y\)に与える影響も閾値周辺で不連続に変化するとは考え難い
  • したがって，この閾値において所得に明確なジャンプがある理由は大卒の効果以外には考えられない

2.2 Rで実践

RDDは閾値周辺のデータを使った回帰分析である局所線形回帰（local linear regression）などを使う。

# 条件2の大卒者と非大卒者のみを取り出す
df_rdd = df %>% filter(is_college2 == 1 | is_college == 0)

# 閾値周辺のサンプルの取り出し
h = 1  # バンド幅：閾値周辺のどの程度の幅のデータを使うか
c = 180 # 閾値
df_local = df_rdd %>% filter(c - h <= score & score < c + h)

取り出した閾値周辺のデータの分布は次のようになっている

# plot
ggplot(df_local, aes(x = score, y = income))+
  geom_point()

# 局所線形回帰
llr <- lm(income ~ score + is_college, data = df_local) 
stargazer::stargazer(llr, type = "text")
## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               income           
## -----------------------------------------------
## score                         13.734           
##                              (12.486)          
##                                                
## is_college                  534.456***         
##                              (14.439)          
##                                                
## Constant                    -1,581.780         
##                             (2,241.268)        
##                                                
## -----------------------------------------------
## Observations                   2,754           
## R2                             0.674           
## Adjusted R2                    0.673           
## Residual Std. Error     190.169 (df = 2751)    
## F Statistic         2,838.943*** (df = 2; 2751)
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

比較的500に近い値を推定することができた